Random-phase-approximation-based correlation energy functionals: Benchmark results for atoms 
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The random phase approximation (RPA) for the correlation energy functional of density functional theory 
has recently attracted renewed interest. Formulated in terms of the Kohn-Sham (KS) orbitals and eigenvalues, it 
promises to resolve some of the fundamental limitations of the local density and generalized gradient approxi- 
mations, as for instance their inability to account for dispersion forces. First results for atoms, however, indicate 
that the RPA overestimates correlation effects as much as the orbital-dependent functional obtained by a second 
order perturbation expansion on the basis of the KS Hamiltonian. In this contribution, three simple extensions 
of the RPA are examined, (a) its augmentation by an LDA for short-range correlation, (b) its combination with 
the second order exchange term, and (c) its combination with a partial resummation of the perturbation series 
including the second order exchange. It is found that the ground state and correlation energies as well as the 
ionization potentials resulting from the extensions (a) and (c) for closed sub-shell atoms are clearly superior to 
those obtained with the unmodified RPA. Quite some effort is made to ensure highly converged RPA data, so 
that the results may serve as benchmark data. The numerical techniques developed in this context, in particular 
for the inherent frequency integration, should also be useful for applications of RPA-type functionals to more 
complex systems. 
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I. INTRODUCTION 

Recent years have seen a revival of interest in the 
random phase approximation (RPA) and its extensions, 
both in the framework of Kohn-Sham density functional 
theory (KS-DFT)'-2-3AS,V,8i9iio,ii,i2,i3,i4,i5,i6,i7 ^j^^in 

Green's function-based many-body theory for ground state 
properties, i** '''^"^! ^! within KS-DFT, the RPA for the en- 
ergy and response function of the homogeneous electron gas 
played an important role in the development of the local- 
density approximation (LDA) as well as the generalized gradi- 
ent approximation (GGA) for the exchange-correlation (XC) 
energy functional E^c^^'^'^ Current interest in the RPA is stim- 
ulated by the concept of orbital-dependent (implicit) XC func- 
tionals, in which the XC energy is represented in terms of 
the KS orbitals and eigenenergies . 25,26,27.28.2 9 wjfjjjjj ^jjjs ^p. 
proach an RPA-type correlation energy functional is most eas- 
ily formulated on the basis of the KS response function. Com- 
pared to LDA/GGA-type explicit XC functionals, implicit 
functionals have several attractive features: (1) the exchange 
can be treated exactly, leading to exchange energies and po- 
tentials which are free of self-interaction;^*' (2) the long-range 
dispersion interaction can be correctly described j'''^'^''^^ (3) 
static coiTelation effects can be incorporated even within a 
spin-unpolarized formalism.'^ 

A systematic formulation of orbital-dependent XC func- 
tionals is possible within KS-based many-body theory, i.e. 
by using the KS Hamiltonian as non-interacting reference 
Hamiltonian in the framework of standard many-body the- 
ory (KS-MBT).Sii2i21 In this approach the exact exchange 
of DFT emerges as first order contribution to a straightfor- 
ward perturbation expansion in powers of e^. All higher or- 
der terms constitute the DFT coiTelation energy. The low- 
est order congelation contribution resulting from perturba- 

(2) 

tion theory, Ec , has been extensively studied for atoms 



and small molecules 



31.32.35.36.37.38 



(2) 

Ec correctly accounts for 



the dispersion interaction^ ''^^ and the corresponding correla- 
tion potential Vc reproduces the shell-structure and asymp- 
totic behavior of atomic coiTelation potentials. ^'^ On the other 
hand, the magnitude of the energies and potentials resulting 
from Ec overestimates the corresponding exact data signif- 

(2) 

icantly. Moreover, Ec is found to be variationally insta- 
ble for systems with a very small energy gap between the 
highest occupied and the lowest unoccupied molecular or- 
bital (HOMO-LUMO gap)"'^^ (as, for instance, the beryl- 
lium atom) and fails to describe chemical bonding in such el- 
ementary molecules as the nitrogen dimer^^ The variational 

(2) 

instability of Ec can be resolved by resummation of suit- 
able higher order contributions to infinite order (e.g. in the 
form of Feynman diagrams). The simplest functional of this 
type is obtained by resummation of selected ladder-type dia- 
grams, i.e. the Epstein-Nesbet(EN)-diagrams. The resulting 
functional is not only found to be variationally stable for all 
neutral and singly-ionized atoms, but also gives more accurate 
correlation energies and potentials than E^^\^'^ However, EN- 
type functionals still face fundamental problems in the case of 
degenerate or near-degenerate systems. A more suitable par- 
tial resummation scheme is needed to establish a universally 
applicable, implicit XC functional, the RPA and its extensions 
being the prime candidates. 

In standard many-body theory, the RPA is obtained by re- 
summation of the so-called ring diagrams.""' This concept can 
be directly transfered to the framework of KS-MBT."^' On the 
other hand, in the context of DFT, the RPA can also be derived 
from the adiabatic connection fluctuation-dissipation (ACFD) 
theorem.^^ The ACFD formalism is, for instance, the con- 
ceptual starting point for the recent development of van-der- 
Waals DFT. ' It has also been applied directly to various sys- 
tems, including jellium surfaces and slabs,-', atoms^' '*-', small 
molecules'*'"'''^'-- and solidsjiiiii All these calculations have 
demonstrated promising features of RPA-based functionals. 
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On the other hand, the results for atoms, ^'•'^^ for which rigor- 
ous benchmark data are available, indicate that the pure RPA 
overestimates correlation energies and potentials as much as 

One is therefore led to consider extensions of the RPA. The 
most obvious starting point for extension is the inclusion of 
the second order exchange (SOX) contribution. However, in 
its pure form it neglects the screeening of the Coulomb in- 
teraction, which is the core feature of the RPA. One would 
thus expect an imbalance between direct and exchange con- 
tributions, when combining the RPA with the pure SOX term. 
A fully screened form of the SOX is easily formulated, fol- 
lowing the line of thought used for the derivation of GGAs.**"* 
The resulting functional, however, is computationally much 
more demanding than the RPA. For that reason it is worth- 
while to examine alternative modifications of the SOX term 
which reduce its net contribution. Given the success of the 

(2) 

EN-resummation in the context of the complete Ec , an EN- 
extension of the SOX term suggests itself (this functional is 
denoted as RSOX in the following). 

The SOX term, be it screened or not, is inherently a short- 
range contribution. This raises the question whether it is suffi- 
cient to account for the complete screened SOX in an approx- 
imate fashion, relying on the LDA. In fact, using this strategy, 
one can easily include all short-range correlation effects be- 
yond the RPA.'*'^ Clearly, the resulting LDA-type functional 
(here labelled as RPA-n) is even more efficient than the RSOX. 

In this work, we study the RPA and these simple extensions 
for a series of prototype atoms and ions, for which highly ac- 
curate reference data are available for comparison. In order to 
provide benchmark results a numerically exact, i.e. basis-set- 
free, approach is used and considerable emphasis is placed on 
all convergence issues involved. As a byproduct of this strive 
for accuracy, a highly efficient scheme for performing the fre- 
quency integration inherent in all RPA-type functionals has 
been developed. This procedure should be useful for applica- 
tions to more complex systems, for which utilizing more than 
the minimum number of grid points for the frequency integra- 
tion would be too demanding. 

A complete implementation of any XC functional requires 
not only the evaluation of the XC energy, but also the in- 
clusion of the corresponding XC potential w^c in the self- 
consistent calculation. The latter step is quite challenging 
in the case of orbital-dependent XC functionals, for which 
t^xc has to be determined indirectly via the Optimized Po- 
tential Method (OPM),2^i2^i21ZMi and, in particular, for RPA- 
type functionals i^^'^^''^^i"^^''*^i^^'"^^i^'' Recently, Hellgren and von 
Barth have reported the first self-consistent RPA correlation 
potentials for spherical atoms, obtained by solution of the lin- 
earized Sham-Schluter equation-- at the GW level4^ How- 
ever, as indicated earlier, the RPA is not consistently improv- 
es) 

ing atomic correlation potentials over Ec ■ In the present 
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work we therefore focus on the perturbative evaluation of all 
RPA-type energies, utilizing self-consistent exchange-only or- 
bitals and eigenvalues. As we will show, the RPA correla- 
tion energy is rather insensitive to the KS orbitals used for its 
evaluation, which clearly supports this perturbative approach. 
This feature, if true in general, will be very important for the 
application of RPA-type functionals to more complicated sys- 
tems, for which a self-consistent implementation is not feasi- 
ble anyway. 

The paper is organized as follows. In the Section II, first the 
RPA correlation energy is formulated in the framework of the 
ACFD formalism. In addition, the general result is reduced to 
an expression valid for spherical systems. In Section III, var- 
ious numerical aspects are discussed, addressing in particular 
questions of accuracy. In Section IV, the RPA correlation en- 
ergies for a number of prototype atoms and ions (we will no 
longer distinguish between neutral atoms and atomic ions in 
the following) are presented and compared to the correspond- 
ing exact data. Section V provides a summary. Atomic units 
are used throughout this paper. 



II. THEORY 

A. RPA correlation energy on basis of the ACFD formalism 

Based on the adiabatic connection and the zero-temperature 
fluctuation-dissipation theorem,"*^'^^ the exact KS correlation 
energy can be written as 

Ec = — — [ du [ d\ I drdr'vcc{r — r') 
27r Jo Jo J 

X [XA(r,r';iu) - Xo(r,r';iu)] , (1) 

where Wcc(r, r') ~ l/|r — r'| is the bare Coulomb interaction, 
Xo is the KS response function, 

Xo(r,r : — : l-c.c., (2) 

^ tu + ei-Sa 

and with As [0, 1], is the density-density response func- 
tion of a fictitious system in which electrons interact with 
a scaled Coulomb potential Auoc(i", i"')' ^nd simultaneously 
move in a modified external potential, chosen such that the 
electron density remains identical to that of the fully interact- 
ing system with A = 1 . Throughout this paper we use the con- 
vention that . . . denote occupied (hole) KS states, while 
a,b, . . . are used for unoccupied (particle) states and p,q, . . . 
for the general case. x\ is related to xo by ^ Dyson-like inte- 
gral equationr^ 



Xx{r,r',iu)^Xo{r,r',iu) + J dri J dr2Xoir, ri, iu)Kx{ri, r2, iu)x\{r2, r\iu) , 



(3) 
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where 

Kx{ri,r2,iu) = Auco(ri, + /4(ri, r2, (4) 

is the Coulomb and XC kernel. 

The RPA correlation energy is obtained from Eq.([T]) if one 
neglects the XC contribution to the right-hand side of Eq.©. 
Integrating over A one ends up with 



E. 



RPA 



du 
2^ 



Tr {Ln (1 - xo{iu)vcc) + Xo{'iu)vcc} 



(5) 

where the trace indicates integration over all spatial coordi- 
nates. It is often more instructive to rewrite the integrand in 
Eq. (|5]l, denoted as Ec{iu), as a power series in the Coulomb 
interaction. 



n=2 



(6) 



B. Correlation energy beyond RPA 

The ACFD theorem provides a natural starting point for the 
development of correlation functionals beyond the RPA: In- 
clusion of some approximation for f^c in the Dyson equation 
(O automatically yields an extension of the RPA. Several ap- 
proximate XC kernels have been introduced in the context of 
time-dependent DFT (TDDFT)i^i^ It is, however, not clear 
whether an approximate f^c designed to provide a good de- 
scription of excited states within TDDFT also leads to accu- 
rate ground state correlation energies. 

A more straightforward extension of the RPA, the so- 
called RPAh- approach, had been proposed by Perdew and 
coworkers.'*'*' They observed that the RPA provides a quite ac- 
curate description of long-range correlation, but is inadequate 
for short-range correlations. On the other hand, the latter can 
be very well approximated by a local or semi-local density- 
based functional (LDA- or GGA-type), 



pRPA+ pLDA 
Ec - ^c,sr 



E, 



RPA 



(7) 



where the LDA for the short-range contribution -Ec,si can be 
obtained by subtraction of the RPA-limit from the full LDA 
correlation energy. 



77.LDA pLDA 77LDA-RPA 

ErST ~ Ec - E^ 



(8) 



This approach is supported by the fact that the gradient cor- 
rection to the short-range correlation is much smaller than that 
to the complete correlation energy."* Though the RPAh- func- 
tional has been used recently to describe the inter-layer disper- 
sion interaction in boron nitride,'^ a direct comparison of the 
RPAh- with exact results is still missing even for closed-shell 
atoms. Using numerically exact RPA correlation energy avail- 
able for atoms, we are able to give a unambiguous assessment 
of the quality of the RPA+ correlation functional. 

In the language of Feynman diagrams, the RPA correlation 
energy is obtained from the second order direct diagram by 



replacing the bare Coulomb interaction by the dynamically 
screened Coulomb interaction. The dominant contribution 
that is missing in the RPA is the second order exchange di- 
agram (SOX), 



E. 



SOX ^ _i (^i 1 1 ai>) {ab \ \ ji) 
2 + £j -Sa-£b 



(9) 



where the notation 



ipq\\rs) ^ I dr I dr -— 



(10) 



has been used for the KS Slater integral. Combining the RPA 
with the SOX term, one obtains a new functional, denoted 
as RPA+SOX. However, one would expect the SOX to over- 
correct the error in the RPA, since the Coulomb interaction en- 
ters the SOX term in its bare, i.e. un-screened, form. Screen- 
ing can be introduced into the SOX term in a systematic way 
by use of the same, dynamically screened interaction as in the 
direct term."*"* Unfortunately, the resulting functional is com- 
putationally much more demanding than the RPA expression. 
A technically much simpler way to reduce the SOX contri- 
bution has been suggested in the context of the second order 
functional Ec'^'' ■.^'^•^^ The inclusion of the direct hole-hole con- 
tribution to the Epstein-Nesbet-type ladder diagrams into the 
SOX term substantially improves second order energies and 
potentials, without introducing any additional computational 
effort. Although the physical background of these ladder di- 
agrams is quite different from dynamical screening, it seems 
worthwhile to analyze this "effective" screening. The result- 
ing correction will be denoted as RSOX, 



E. 



RSOX 



2 ^ 

i,j,a,b 



{ij\\ab) iab\\jt) 



£b - {ij 1 1 ij) 



(11) 



C. RPA correlation functional for spherical systems 

In the case of spherical systems, the KS potential only de- 
pends on the radial coordinate r, Vs{r) — Vs{r), and each KS 
orbital can be written as the product of a radial orbital and a 
spherical harmonic Y/m(0, ip), 

M^) ^ <Pnimir) = ^^Yi^ie,^), , (12) 
r 

where n, I and m are the principle, angular and magnetic 
quantum numbers, respectively. The Pni are solutions of the 
radial KS equation. 



1 / d^ l{l + l) 

2 V dr^ 



Pnl{r) = EnlPnlir) 



(13) 

Two-body functions like the Coulomb interaction and xo 
can also be decomposed according to the spherical symme- 
try. To simplify notations, we use the following decomposi- 
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tion convention. The Coulomb interaction is expanded as 



" 

L 

^ rLM(0,^)>7A/(^'>^') , (14) 



L=0 



where 



VL{r,r') ■.= r^/r!^+^ (15) 



with r< = Min(r, r') and r> = Max(r, r'). The response 
function xo can be written as 



, ^ 2i + 1 X0L{r,r',iu) 



L=0 



A/=-L 



The L-dependent radial response function xol(^j can 
be calculated utilizing the radial orbitals 

XOL{r,r',iu) = -YCL-iaaD^aa{u)^zaa{r)^iaa{r') , 
iacr 

(17) 



where 



Diaa{u) 



(2Z, + l)(2Za + l) fh la L 



2L + 1 V 



--acr ^1(7 1 



(18) 

(19) 
(20) 



Using the multipole expansion of both v^e and xo^ the 
building block of the RPA correlation energy, Tr{xow}, is 
be obtained as 

POO 

Tr {Xov} = ^(2i + 1) / drdr'xoLir, r' , iu)vL{r' ,r) . 
L 

(21) 

There are two options for the calculation of the radial integral 
in Eq. (IHT i: 

a. Real space approach: In this approach, Xol(?', is 
calculated on a discrete radial mesh, which allows to evaluate 
(I2TI 1 by direct numerical integration, 

Tr {xqv} 

= Yi'^L + l)Yw{n)w{rj)xoL{n,rj)vL{rj,n) 



= ^(2L + l)Tr{xoL«L} , 



(22) 



where Wi denotes the radial integral weight at mesh point i. In 
case of powers of Tr {xo"*^} one has 

Tr{(xo«n=5](2L + l)Tr{(xoL«in . (23) 

L 

The sum over n in Eq. ^ then leads to 

ESu)^Y.^2L + l)Tr{LYv{l-xoLVL)+XoLVL] ■ 

(24) 

b. Orbital-product space approach: In this second ap- 
proach one inserts (fT6] l and ( fTTI l into Eq. (ISTT i. Using the 
radial Slater integrals 

RL;iaa,3ba' dr dr' <^iaa{r)v L{r , r')^ jba' {r' ) 



"'0 



one finds 

Tr{xow} = y^(2L + 1) 



Jba' {■!■ ) , 
(25) 



^ J dr J dr'<i>,aa{r)vL{r,r')^,a,y{r') 

- ^(2L + 1) XI CL-taaD^aa{u)RL-ra 



(26) 



With the definitions 



^;iaCTj&(T' C [^■iaaRL;iaa,jb(j' \/ CL-jba' (27) 
SL:iaa,jba' ■= — \/ Diaa {u)VL-iaa,jba' \J Djb„'{u) (28) 

one ends up with 

Tr{xoi^} = Y.^2L+1)Y,Sl 

L iaa 

= X(2L + l)Tr{5L} . (29) 

L 

One can furthermore show that 

Tr {(xo^^)"} = + l)Tr {(^l)"} (30) 



so that 



£;,M = -X(2L + l)TrJX 



, Ji=2 



= X(2L + l)Tr{Ln(l-5L) + 5L} 
i 

= X(2L + l)[Ln (Dot (1 - Sl)) + Tr {5l}](31) 

L 

The final expressions for Ef^^ are quite similar in the real- 
space and orbital-product-space approaches, but their numer- 
ical efficiency can be very different, depending on the size of 
the system. In the real-space approach, the dimension of the 
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matrix involved is determined by the number of mesh points 
used for radial integration, /max- As /max is never larger than 
a few thousand even for heavy atoms, the resulting memory 
requirement is quite low. On the other hand, xol needs to 
be constructed on the radial mesh for each frequency, which 
can be very cpu- time-intensive. The situation is quite different 
in the case of the orbital-product space approach. Here the di- 
mension of the matrix involved is given by iVocc x ^vir, where 
A'occ denotes the number of occupied orbitals and A'vir is the 
number of unoccupied orbitals taken into account. A'occ x ^vir 
can be easily as high as tens of thousands. However, the ma- 
trix Vl, Eq.(l27]l. is independent of frequency and can be cal- 
culated in advance and stored in the memory. Limitations of 
the available memory can be circumvented by taking advan- 
tage of the fact that, for given L, Vl and Sl are quite sparse 
(usually the ratio of non-zero elements is less than 1%). One 
can therefore use standard sparse matrix techniques to reduce 
the storage requirement and accelerate the computation of the 
determinant. 



III. NUMERICAL DETAILS 
A. Hard-wall cavity approach 

The RPA correlation energy depends on the occupied as 
well as on the unoccupied KS states. For free atoms, the 
spectrum of the unoccupied states includes both discrete Ry- 
dberg states and continuum states. However, the handling of 
continuum states in the evaluation of the correlation energy 
is nontrivial.^^. Moreover, the presence of continuum states 
causes additional problems in the context of orbital-dependent 
functionals: One does no longer find a solution of the corre- 
sponding 0PM equation which satisfies the standard bound- 
ary conditions for the correlation potential. To resolve this 
problem, we developed a hard-wall cavity approach, ^^"^^i^ in 
which the KS equation is solved on a discrete radial mesh with 
hard-wall boundary conditions imposed at a finite, but large 
radius Rmax- The same approach is used in the present work. 

Its crucial parameters are the cavity radius /?max as well as 
the energetically highest state (characterized by its principle 
quantum number nmax) and the highest angular momentum 
^max included in sums over virtual states. In the following, 
neutral Ar is used for a systematic convergence study with 

respect to /?max, "max and ^max- 

We first consider /?max- Rmax has to be chosen so large, 
that all ground state properties, and, in particular, the corre- 
lation energy, do no longer change when /?max is increased 
further. However, any increase of /?max directly affects the 
spectrum of the positive energy states, i.e. the density of states. 
In order to keep the space available for virtual excitations con- 
stant, when increasing Rmax, one therefore has to fix the en- 
ergy Emax of the highest unoccupied state rimax taken into ac- 
count. In the case of very high-lying virtual states one has a 
simple relation between Rmax and Emax, resulting from the 
fact that high-lying states are no longer sensitive to the de- 



TABLE I: Convergence of , the exact exchange energy and the 
eigenvalue of the highest occupied KS orbital (ehomo) obtained by 
self-consistent exchange-only calculations for Ar for different cavity 
radii i?max (with nmax/-Rmax = 10 Bohr"^ and /max = 4, i?max in 
Bohr, all energies in Hartree). 



Rn 



RPA 



ehomo 



5 1.0023 30.2059 0.5772 
8 1.0027 30.1749 0.5909 
10 1.0028 30.1747 0.5908 



TABLE II: Convergence of full E^^^ (Column 3) and E^'^^ within 
the frozen core approximation excluding virtual excitations of the 
Is, 2s and 2p electrons (Column 4) of Ar with respect to rimax (with 
fimax = 10 Bohr and /max = 4, all energies in Hartree). 



^max 


£max 




-E^^\FC) 


25 


25.1 


0.6840 


0.3961 


50 


111.9 


0.9097 


0.3980 


100 


471.2 


1.0028 


0.3980 


150 


1077.6 


1.0273 


0.3980 


200 


1930.8 


1.0354 


0.3980 


250 


3030.9 


1.0384 


0.3980 


300 


4377.6 


1.0398 


0.3980 


350 


5951.9 


1.0399 




400 


7754.8 


1.0400 





tailed structure of Vs, 



Rn 



(32) 



The space available for virtual excitations is therefore kept 
constant, as soon as the ratio n-maxl Rmax is fixed. Table 
11] shows the values of E?~^^ for Ar obtained with different 



TABLE III: Convergence of full B^^* (Column 3) and Ef^'^ within 
the frozen core approximation excluding virtual excitations of the 
Is, 2s and 2p electrons (Column 4) of Ar with respect to /max (with 
Rmax = 10 Bohr and Wmax = 100, all energies in Hartree). 



^max 


^RPA 


-£;?PA(FC) 


2 


0.7661 


0.2875 


4 


1.0028 


0.3980 


6 


1.0431 


0.4185 


8 


1.0529 


0.4241 


10 


1.0562 


0.4262 


12 


1.0574 


0.4271 


14 


1.0580 




16 


1.0582 
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^max, but fixed nmax/i?n 



10 Bohr , which corre- 



sponds to an energy cut-off of about 500 Hartree. For com- 
parison, the corresponding exchange energy and the eigen- 
value of the highest occupied orbital, enoMO^ resulting from 
exchange-only calculations are also listed. One observes that 
jTjRPA jggg sensitive to i?max than the exchange energy, 
which is consistent with the fact that the length scale related 
to the RPA correlation energy is smaller compared to that of 
the exchange. 

Argon is the heaviest atom considered in this work. We 
have also made systematic convergence tests for other, less 
compact atoms like Na and Mg. For all atoms considered in 
this work, the choice i?max = 10 Bohr leads to errors less 
than 1 mHartree. 

With i?,nax fixed, one can now examine the convergence 
of E^^^ with respect to Umax and /max- Tables and |lll] 
show Ef-^^ for Ar obtained with different rimax and Z,nax- In 
general, the absolute value of Ef^^ converges quite slowly 
with respect to both parameters. The slow convergence with 
respect to n,„ax mainly originates from the innermost shell 
— unoccupied states with high energies are only important 
for the description of virtual excitations of the highly local- 
ized, low-lying core states. In practice, fortunately only en- 
ergy differences related to the valence electrons are really 
relevant. One would thus expect to achieve high accuracy 
for these energy differences with much more moderate val- 
ues for rimax-This suggests to rely on the frozen core (FC) 
approximation, in which excitations from core levels are ex- 
cluded. Tables HIl and Hill demonstrate that the FC approxima- 
tion for Ef^^ converges much faster with increasing rimax- 
Even for a quite moderate n„iax of 25, corresponding to 
£max ~ 25 Hartree, Ef-^^ is already converged to an accu- 
racy of 2 mHartree. On the other hand, the convergence be- 
havior of Ef^^ with respect to /,„ax is not improved by the 
FC approximation. As one of the main aims of this work is to 
provide benchmark results for a set of prototype atoms, most 
results reported in this work are obtained without evoking the 
FC approximation. The results reported in the next section are 
obtained for n„iax = 300 and /max = 14, which ensures an 
accuracy of 1 mHartree for Ar and better for all lighter atoms. 



B. Frequency integration 

Any calculation of RPA energies involves two time- 
consuming steps: The first is the evaluation of all Slater in- 
tegrals involved, i.e. of the matrix Rl, Eq.(l25]l. In the present 
work the Slater integrals are calculated by numerical integra- 
tion on the radial grid, using standard finite differences meth- 
ods. Once Rl is available, the second step is performing the 
frequency integration in Eq.(|5]l. In order to understand the 
most appropriate way to do this frequency integration let us 
consider the integrand for some prototype atoms. Figure [T] 
shows u^Ec{iu) for He, demonstrating the fact that Ec{iu) 
falls off as u^^ for extremely large u. This behavior can 
be easily understood on the basis of Eq. (fT9] i: For frequen- 
cies beyond the maximum excitation energy Se — e™^'^ — 
included in the calculation (or provided by the basis set) 
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FIG. 1: (Color online) u^Ec{iu) vs u for He (maximum excitation 
energy Se ~ 4000 Hartree). 
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FIG. 2: (Color online) u^Ec(iu) and u'^Ec 
large u for the case of He. 



vs u for moderately 



Diaa{u) and thus Sl{u) decay as which allows a pertur- 
bative evaluation of (|3TI ) in powers of Sl (u), with the second 
order term dominating the resulting Ec{iu). 

On the other hand, for the more important range of large 
frequencies less than 5e a decay close to u^^ is found, as 
shown in Figure |2l The same behavior is observed for each 
individual shell, as illustrated by the Ec{iu) obtained by ex- 
citation of only the A/-shell of neutral Ar, also included in 
Figure |2] 

This power law decay suggests a transformation of the fre- 
quency interval < u < oo to some finite interval via a power 
law transformation, as for instance 



1 



1 + (u/s) 



< a; < 1 , 



(33) 



7 




X 



FIG. 3: (Color online) Integrand Ec{iu)du/dx of frequency integral 
after the transformation i33i vs x for He, Ar^^^ as well as the M- 
shell of neutral Ar. 

giving more weight to large u than an exponential transforma- 
tion. The scale factor s is intrinsically related to the minimum 
energy required for a virtual excitation, which is roughly given 
by the eigenvalue difference elumo — (e) ~ for a shell 
with average eigenvalue (e) . The most appropriate s can only 
be determined empirically. For all atoms considered in detail 
the choice s — 2|(e)| seemed to work reasonably well (see 
also below). The result of the transformation ( [33] l is shown 
in Figure [3] for He, Ar^^+ as well as the A/-shell of neutral 
Ar One obtains a smooth function of x, with values remain- 
ing on the same order of magnitude for all x. This ensures 
a rapid convergence of the numerical integration over x with 
the number of grid points. 

However, the frequency integration in Q suffers from the 
fact that each shell in an atom (or molecule) introduces a new 
energy scale for the virtual excitations. This is most easily 
verified by plotting uEc{iu) on a double logarithmic scale, as 
done in Figure|4]for He, Ne and Ar. Figure|4]demonstrates that 
there are two relevant scales for Ne, three in the case of Ar In 
fact, the plots confirm that the behavior of Ec{iu) changes at 
roughly twice the average eigenvalues of the shells involved, 
with the position of the highest energy transition point being 
somewhat less pronounced (ex(Ne)=— 30.8, eL(Ne)=— 1.1; 
ej^(Ar)=-114.4, ei,(Ar)=-9.4, eji/(Ar)=~0.7 — all values 
in Hartree). It is therefore necessary to split the frequency 
integration from to oo into intervals associated with these 
individual energy scales. Let us call the boundaries of the in- 
tervals bi, 

= bo < bi < . . . < bn = oo (34) 

with n denoting the number of shells. The intervals are cho- 
sen such that the characteristic energy scale s„ of shell n is 
bracketed, 

b^-l < s, < 6, . (35) 




u [Hartree] 

FIG. 4: (Color online) uEdiu) vs u for Ne and Ar. 



In practice, Si = 2|(ei)| and bi = 4|(ei)| seemed to provide a 
reasonable partioning of the complete frequency range. Eq.(|5| 
may then be decomposed as 

/ Ec{iu)du = Ec{iu)du . (36) 

In order to account for the piecewise decay of Ec{iu) the 
transformation 

^ _ [1 + Si/{bi - u - bj^i ^^^^ 

[1 + (u - 6j_i)/si] 
^ ^ [l + s,/{b, ^^g^ 
dxi ' [1 + Si/(6i - - 

{i — 1, . . .n) is most suitable. Eq.(l36b can then be written as 

/ Ec{iu)du = / dxi-—Ec{iu{xi)) . (39) 
Jo 7^1 J "-^^ 

The success of this frequency partioning plus power law trans- 
formation scheme is demonstrated in Figure |5] in which the 
final integrands du/dxi Ec{iu{xi)) are plotted for neutral 
Ar In all three energy regimes a rather smooth integrand 
is obtained, which allows the application of Gauss-Legendre 
quadrature to each interval. As a result, the error obtained for 
a given total number of grid points Nu used for the Gauss- 
Legendre quadrature is rather small already for very low Nu, 
as illustrated in Figure|6] The most critical interval in Eq.([39l) 
is the highest energy range, covering in particular excitations 
of the Is-state. For that reason the eiTor is even lower if only 
excitations of the valence shell are included (i.e. in the FC ap- 
proximation), as is obvious from the error found for He or the 
M -shell of neutral Ar 
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FIG. 5: (Color online) Integrand Ec(iu(x))du/ dx of partioned fre- 
quency integral ( 139b resulting from the transformation ( I37t for neu- 
tral Ar. 




FIG. 6: (Color online) Absolute error resulting from Eas.(l39t and 
J37b as a function the total number of grid points Nu used for the 
Gauss-Legendre quadrature (note: A'^„ is the sum of the number of 
grid points used in the individual intervals). 



IV. RESULTS 

A. Sensitivity to form of KS orbitals 

Standard KS-DFT calculations are based on the self- 
consistent solution of the KS equations, which requires the 
evaluation of the XC potential Wxc(r) = 5E^c[n]/ 5n{Y). In 
the case of LDA and GGA functionals the calculation of 
?;xc (r) is straightforward. On the other hand, a self-consistent 
implementation represents a much more serious problem for 
RPA-type functionals. First of all, in the case of orbital- 
dependent functionals Uxc has to be determined via the OPM, 



i.e. by solution of an integral equation. The solution of the 
OPM integral equation is well-established for the exact ex- 
change and managable, though rather intricate, for the second 
order correlation functional iJ^-'j^^iSi^^^ Its implementation 
for RPA-type functionals, however, is much more challenging, 
so that a full solution has only been reported very recently,— 
On the other hand, a self-consistent implementation is only 
then advantageous if the resulting correlation potential leads 
to an improved total KS potential. It has been demonstrated 
that this is not the case for some standard GGAs''" and for 
^(2)^38^ which by far overestimates correlation effects. The 
situation is not yet completely clear in the case of the RPA. 
The first RPA potentials available'*^' seem to improve over 
Vc — SEc /Sn{r) in the asymptotic regime, but otherwise 

(2) 

often follow Vc ■ 

Independently of this more fundamental aspect, one might 
ask whether a self-consistent implementation is really neces- 
sary to obtain accurate RPA correlation and thus ground state 
energies. Clearly, a purely perturbative treatment of RPA- 
type functionals on the basis of a self-consistent calculation 
with the exact exchange would allow their application to much 
more complex systems, for which the solution of the coiTe- 
sponding OPM integral equation is beyond current computer 
resources. In fact, experience with conventional density func- 
tionals suggests that, at least for atomic systems, the RPA 
correlation energy is not sensitive to the detailed structure of 
t'xc(r). 

In order to verify this expectation, the RPA ground state 
energy (i.e. the sum of the KS kinetic energy, the Hartree 
term, the exact exchange and E^^^) has been calculated by 
insertion of KS orbitals resulting from different XC function- 
als: Orbitals obtained from self-consistent calculations with 
only the exact exchange (EXX-only), but neglecting Vc com- 
pletely, are compared with self-consistent LDA and GGA 
orbitals. The results for a number of atoms are collected 
in Table |IV] which also includes recent self-consistent RPA 
energies,'*^ whenever available. Table II V I confirms the expec- 
tation: Even though the KS potentials obtained by EXX-only 
calculations differ substantially from their LDA and GGA 
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TABLE IV: Absolute RPA total energies (in Hartree) resulting from insertion of different KS orbitals. The last row lists the self-consistent 
RPA total energies given in Ref 43. 



KS orbitals 


He 


Be 


Ne 


Mg 


Ar 


N 


Na 


LDA 


2.945 


14.751 


129.140 


200.293 


527.905 


54.735 


162.475 


BLYP 


2.944 


14.752 


129.142 


200.296 


527.910 


54.737 


162.478 


EXX-only 


2.945 


14.752 


129.143 


200.298 


527.913 


54.738 


162.480 


RPA 


2.945 


14.754 


129.143 


200.296 


527.908 







counterparts, the differences between the resulting RPA en- 
ergies are quite small. The same is true for the deviations be- 
tween the perturbative RPA energies on EXX-only basis and 
the self-consistent RPA results. This result is expected to hold 
quite generally, as long as one does not examine a quantity 
which is particularly sensitive to the correlation potential (as, 
for instance, the quantum defect of high Rydberg states''' ). In 
fact. Table |IV] indicates that even a perturbative treatment of 
both the exact exchange and the RPA may be legitimate for 
very complex systems, in which even self-consistent calcula- 
tions with the exact exchange are too expensive. In the fol- 
lowing sections, self-consistent EXX-only orbitals are always 
used for the evaluation of the RPA correlation energy. 



B. RPA correlation energies of spherical atoms 

In this work, we focus on atoms with closed or half-filled 
shells, for which the ground state KS potentials for the two 
spin-channels are both spherically symmetric. Table [Vl lists 
the correlation energies obtained from all four RPA-based 
functionals for a series of atoms. To see trends more clearly, 
the relative errors resulting from the various functionals with 
respect to the exact correlation energies^^ are plotted in Fig- 
ure [T] Not surprisingly, the pure RPA always overestimates 
the true correlation energy. Adding the short-range correc- 
tion within the LDA (RPA-n) improves the results remarkably, 
reducing the mean absolute error by more than an order of 
magnitude. As expected, the unscreened SOX contribution 
by far overcorrects the error of the RPA. On the other hand, 
the inclusion of EN-corrections into the SOX term (RSOX) 
reduces this overcorrection significantly. In general, both the 
RPAh- and the RPAh-RSOX produce more accurate correlation 
energies than the LYP-GGA, at least for the set of atoms con- 
sidered in this work. Moreover, for light atoms one observes 
a tendency of the RPAh-RSOX to be superior to the RPA-n. 



C. Ionization potentials 

Even more important than the accuracy of total (correlation) 
energies is the accuracy of energy differences. In the case of 
atoms the ionization potential (IP) serves as the prototype en- 
ergy difference for assessing the quality of any approximation. 
Much more than the total atomic Ec, the IP probes the descrip- 
tion of the correlation of the valence states. We have therefore 
calculated the IPs resulting from the four RPA-based correla- 



TABLE V: Absolute correlation energies (in Hartree) of closed sub- 
shell atoms calculated from the RPA and RPA+ functionals by in- 
sertion of the exact EXX-only KS orbitals in comparison with exact 
data.^" The last row provides the mean absolute error (MAE) with 
respect to the exact energies. 



atom 


exact 


RPA 


RPA+ 


RPA+SOX RPA+RSOX 


LYP 


He 


0.042 


0.083 


0.047 


0.035 


0.044 


0.044 


Li+ 


0.043 


0.087 


0.048 


0.039 


0.045 


0.048 


Be2+ 


0.044 


0.088 


0.048 


0.041 


0.045 


0.049 


Li 


0.045 


0.112 


0.059 


0.029 


0.053 


0.053 


Be+ 


0.047 


0.122 


0.066 


0.031 


0.056 


0.061 


g2+ 


0.049 


0.131 


0.073 


0.030 


0.059 


0.067 


Be 


0.094 


0.179 


0.108 


0.058 


0.097 


0.095 


B+ 


0.111 


0.205 


0.131 


0.066 


0.110 


0.107 




0.126 


0.228 


0.151 


0.073 


0.123 


0.114 


N 


0.188 


0.335 


0.201 


0.146 


0.178 


0.192 


0+ 


0.194 


0.345 


0.208 


0.155 


0.184 


0.207 


p2+ 


0.199 


0.355 


0.215 


0.162 


0.189 


0.218 


Ne 


0.390 


0.597 


0.400 


0.340 


0.367 


0.384 


Na+ 


0.389 


0.599 


0.398 


0.350 


0.371 


0.399 


Mg2+ 


0.390 


0.601 


0.398 


0.358 


0.375 


0.411 


Na 


0.396 


0.626 


0.410 


0.349 


0.383 


0.408 


Mg+ 


0.400 


0.634 


0.415 


0.360 


0.389 


0.427 


Af+ 


0.405 


0.642 


0.420 


0.369 


0.395 


0.442 


Mg 


0.438 


0.687 


0.453 


0.387 


0.427 


0.459 


A1+ 


0.452 


0.706 


0.468 


0.402 


0.440 


0.481 


Si^+ 


0.463 


0.722 


0.481 


0.414 


0.451 


0.497 


P 


0.540 


0.850 


0.554 


0.482 


0.523 


0.566 


S+ 


0.556 


0.873 


0.573 


0.499 


0.538 


0.588 


Cl2+ 


0.570 


0.893 


0.589 


0.514 


0.552 


0.605 


Ar 


0.722 


1.101 


0.742 


0.661 


0.701 


0.751 


K+ 


0.739 


1.126 


0.763 


0.681 


0.720 


0.771 


Ca^+ 


0.754 


1.150 


0.783 


0.702 


0.737 


0.788 


MAE 




0.196 


0.015 


0.039 


0.011 


0.018 



tion functionals for a number of atoms. In order to avoid any 
uncertainty associated with spherical averaging, the compari- 
son is restricted to atoms for which both the KS potential of 
the neutral ground state and that corresponding to the ionic 
state are spherical. The results are collected in Table I VII The 
most noteworthy features of these data are: (1) The pure RPA, 
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TABLE VI: First ionization potentials (in Hartree) of spherical atoms calculated from total energy differences (IP= iStot {N — 1)) — Etot (N)), 
using different XC functionals. The last row provides the mean absolute error (MAE) with respect to the exact results.— Self-consistent EXX- 
only KS orbitals are used as input orbitals. 



atom 


exact 


EXX 


RPA 


RPA+ 


RPA+SOX RPA+RSOX 


BLYP 


Li 


0.198 


0.195 


0.220 


0.205 


0.185 


0.203 


0.201 


Be+ 


0.669 


0.666 


0.700 


0.683 


0.656 


0.677 


0.681 


Be 


0.343 


0.295 


0.352 


0.338 


0.323 


0.336 


0.329 


B+ 


0.924 


0.861 


0.935 


0.919 


0.897 


0.912 


0.904 


Na 


0.189 


0.179 


0.206 


0.191 


0.178 


0.191 


0.194 


Mg+ 


0.552 


0.540 


0.573 


0.557 


0.543 


0.555 


0.566 


Mg 


0.281 


0.242 


0.294 


0.279 


0.268 


0.279 


0.279 


A1+ 


0.691 


0.643 


0.706 


0.691 


0.677 


0.688 


0.691 


MAE 




0.028 


0.017 


0.005 


0.015 


0.005 


0.009 



though showing significant improvement over the EXX-only 
approximation, generally overestimates the true IPs; (2) Cor- 
respondingly, the RPAh-SOX underestimates IPs (consistent 
with the unscreened nature of the pure SOX term); (3) Both 
the RPAh- and the RPAh-RSOX significantly improve over the 
pure RPA results, and are even more accurate than BLYP, the 
"most accurate" standard GGA. 



V. SUMMARY AND CONCLUSIONS 

In this work, we provide benchmark results for the RPA 
and three simple extensions, allowing for an unambiguous as- 
sessment of the functionals by comparison with exact data. 
Our results confirm earlier observations of the limited appli- 
cability of the pure RPA: The RPA substantially overestimates 
correlation energies, which then results in an overestimation 
of energy differences like ionization potentials. On the other 
hand, our results also demonstrate that already quite simple 
extensions of the RPA can be superior to standard GGAs: 
Adding either short-range corrections within the LDA (RPA-n) 
or a suitably 'screened' second order exchange contribution 
(RPAh-RSOX) significantly improves both absolute energies 
and energy differences. This success is consistent with the ex- 
pectation that the dominant source of error in the RPA func- 
tional is the missing short-range SOX contribution. 

It seems worthwhile emphasizing that the first of these ex- 
tensions, the RPA-H, essentially comes at no cost: Compared 
to the computational demands of an RPA-calculation, the cost 
of the LDA for the non-RPA correlation is irrelevant. More- 
over, a systematic improvement of the RPAh- by inclusion of 



gradient corrections for the non-RPA correlation contributions 
suggests itself. The RPA+RSOX involves an evaluation of the 
orbital-dependent SOX term, which is computationally almost 
as demanding as the calculation of the RPA energy itself, but 
still much less expensive than that of the fully screened SOX 
term. 

In this work also several technical aspects of RPA- 
calculations have been studied systematically, of which two 
should be relevant beyond the regime of atoms considered 
here. The first of these aspects is the sensitivity of the RPA- 
expression to the orbitals and eigenvalues used for its eval- 
uation. It turned out that the character of the KS spectrum 
inserted into the RPA has little impact on the resulting energy. 
Even if KS states obtained by LDA calculations are used the 
deviations from more accurate data remain small. In order 
to cover systems with more than one occupied shell, we have 
developed a partioning scheme for the frequency integration 
inherent in all RPA-type functionals, which, together with a 
suitable transformation of the integration variable, allows to 
perform the frequency integration with a minimum number of 
mesh points. Both these numerical techniques should be par- 
ticularly helpful for applications to more complex systems. 
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